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Abstract 

We study the stability of the (87) Sylvia system and of the neighborhood of its 
two satellites. We use numerical integrations considering the non-sphericity of 
Sylvia, as well as the mutual perturbation of the satellites and the solar pertur- 
bation. Two numerical models have been used, which describe respectively the 
short and long-term evolution of the system. We show that the actual system is 
in a deeply stable zone, but surrounded by both fast and secular chaotic regions 
due to resonances. We then investigate how tidal and BYORP effects modify the 
location of the system over time with respect to the instability zones. Finally, 
we briefly generalize this study to other known triple systems and to satellites 
of asteroids in general, and discuss about their distance from mean-motion and 
evection resonances. 

Keywords: Celestial mechanics. Satellites of asteroids. Resonances, orbital. 
Satellites, dynamics 



1. Introduction 

A large number of satellites of asteroids have been discovered si nce the dis- 



cover y of the satellite Dactyl, thanks to the Galileo flyby of (243) Ida (jBelton et al 
|l996). As of today, there is 206 known systems (binary, triple and quintuple), 
following the Johnston's archive onl ine d atabasjll (see also th e online database! 



described bv lPravec fc Harri^ ( 2007 ) and Pravec et al. ( 2011 )') and it is believed 



that small binaries could r epresent a fraction of 15% of the NEA population 
iMargot et al.. 200l IPravec et al...2006.) . Triple systems are rare and only nine 



known systems have been reported up to now in the entire Solar System. 
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The dynamical evolution and formation mechanisms of these systems are 
highly dependent on the size ratio between the secondaries and the primary. 
If this ratio is very small, as in the case of the Main-Belt asteroid (243) Ida, 
the systems are similar to the classical dynamical p r oblem of a massless satel- 



lite orbiting a planet (see for example iKozai 



19591 1963). this one being re- 



Chauvineau et al ]|l993HScheer"esl 



piacea by a possibly nigniy elongatea eiiipsoia {\^m 

1994t IScheeres et al. 19961 Compere et al. 2012a ). On the other hand, systems 
with similar size components, as the Near-Earth asteroid (66391) 1999 KW4, 
have to be described taking into account both their shapes and their rota- 
tions. A lot of studies have been realized on the expression of the full two- 
body problem and the study of its characteristic s (*Maci eiews kil ll995l: [ScheeresI 
2OO2I: IPahnestock fc ScheeresI I2OO8I: iBoue fc Lask ar 200^^. Similarly, empha- 



sis has been given during the past decade on the description of dissipative 
effects on binary systems, like tidal effects ( Mathis &: Le Poncin-Lafittd 2009 ■ 
Goldreich fc Sarill2009l: pTavlor fc Margotll2010ll201ll) or B YOR P (ICuk fc Burm I 
2005HCuk &: Nesvornvl2010HMcMahon fc Scheeresl2010HSteinberg fc Sarj20rili 



We studied in this paper the dynamics and stabili ty of the systeni (87) S ylvia, 
which was the first triple asteroid system discovered ( Marchis et al.|[2005al ). The 
specificities of this system place it in the first class described above. Sylvia, 
discovered in 1866, is a low-eccentric and midly-inclined asteroid located in the 
outer Main Belt . Its long-term evolution has been investigated th rough the 
AstDys project ( Milani &: Knezevic 1998t lKnezevi6 fc Milanil[2003h giving its 
proper orbital elements (a = 3.486 AU, e = 0.0537, i = 9.85°) and its secular 
fondamental frequencies (n = 55.2977yr, g ^ 134.798" /yr, s = -130.782" /yr). 
Its orbit has been found to be slightly chaotic, exhibiting a Lyapunov time of 
~ 1.4 Myr. 

The two satellites of Sylvia present near-circular and near-equatorial orbits, 
and have a mass ratio of about 10~^ and 10~^ with Sylvia. The outermost 
satellite, Romulus, is approximately ten times more massiv e than the inn e rmost 
one, Remus, for a semi-major axis twice as important. IWinter et al. (120091 ) 
studied the system and found that the satellites could be highly unstable when 
the oblateness of Sylvia (even a small fraction) is not taken into account. Indeed, 
the oblateness of the asteroid, as well as the short distance of the satellites 
from its surface (~ 5 and 10 radius of Sylvia), critically increase the precession 
frequencies and prevent them from commensurabilities with frequencies arising 
from other gravitational perturbations. 

Our aim is the understanding of the dynamical mechanisms present in the 
system and in its neighborhood. We then generalize some of the results to the 
other triple systems, and, in a general way, to the systems similar to (87) Sylvia, 
e.g. with a small size ratio and a primary diameter of the order of ^ 100 km. 
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2. Study of the (87) Sylvia system 



The gravita tional potent ial of Sylvia is modeled by a spherical harmonics 
expansion (e.g. iKaulal . 119661 ): 



t/(r,A,, 



EE 

n— 711—0 



Re 



P„,m(sin (f)) Cn.m COs(toA) + Sn 



sh\{mX) 
_ (1) 

where /i is the gravitational constant of the central body, Re is its radius, (r, A, 0) 
are the spherical coordinates of the satellite, Cn.m and Sn,m are physical con- 
stants depending on the shape of the main body and named coefficients of the 
expansion (n is the degree and m is the order of the coefficients), and Pn,m are 
the associated Legendre polynomials. 

The coefficients of this expansion are computed using the freely available 
software archive SHTOOL^ developed by Mark W i eczorek and using th e con- 
vex shape model of Sylvia (jKaasalainen et al.l 120021 : iMarchis et al.ll2006 ) avail- 
able on DAMn0 (Database of Asteroid Models from Inversion Techniques). 
The asteroid shape models are represented by polyhedrons with triangular sur- 
face facets. Using this shape model as an input of a home-made code (which call 
the functions SHExpandLSQ, MakeGridDH and CilmPlus of SHTOOLS), we 
computed the spherical harmonics coefficients of Sylvia up to the tenth degree 
and order. After having run test integrations, we conclude that a A*^ order and 
degree harmonics development is sufficient to precisely approximate the pertur- 
bation due to the shape of Sylvia. The coefficients up to that degree and order 
are pres ented in Table HI More details on how they have been computed can be 
found in ICompere et al. [ (l2012bh . 

In Table [2j the orbital elements and some parameters used for the integra- 
tions are presented. The incertitudes (when known) are also give n. The orbital 
elements and diameters of the satellites have been taken from iMarchis et al 



( 2005a( ) and their mas ses from Winter et ahl ( 20091 ). The orbital elements of 
the two satellites from Marchis et al. ( 2005al) correspond to different epochs of 
reference, so we took as a reference epoch the mid-time between these epochs 
(JD 2453248) and move the mean anomalies of the satellites to this date, by 
considering fixed mean motions. The heliocentric orbital elements of Sylvia at 
this epoch, as well as its radius and rotation period , have been taken from the 
JPL service. The mass of Sylvia was o btained from Marchis et al.l (l2005al) and 
the ecliptic coordinates of its pole from Drummond fc Christou ( 2008[) . 

Based on this shape model, we used two dynamical models to investigate 
the dynamics of the system. One was used for short-term and precise numerical 
integrations, while the second one was used to study the secular evolution of 
the system. 



''Available at http : //www . ipgp . f r/- wieczor/SHTOOLS/SHTOOLS ■ html | 

*This database is available at http: //a stro .troja.mff ■ cuni ■ cz/projects/asterolds3D/web.php| 
(see iDurech et al.ii2010i for more details) 
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Table 1: Main spherical harmonics coefficients of Sylvia. 



Sylvia 

radius R 130.47 km (c7=6.65) 

mass M 1.478 x lO^^ kg 

rotation period 5.184 li 
ecliptic coordinates of the pole a=100 ° , ^=62 ° 
Remus 

semi-major axis 706 ± 5 km 

eccentricity 0.016 ± 0.011 

inclination 2 ° ± 1 ° 

mean anomaly 96.087 ° 

argument of pericenter 314 ° 

longitude of node 97 ° 

mass 2.154 x 10" kg 

diameter 7 ± 2 km 
Romulus 

semi-major axis 1356 ± 5 km 

eccentricity 0.001 ± 0.001 

inclination 1.7 ° ± 1 ° 

mean anomaly 324.308 ° 

argument of pericenter 273 ° 

longitude of node 101 ° 

mass 3.6625 x lO^^ kg 

diameter 18 ± 4 km 



Table 2: Orbital elements, physical parameters and corresponding incertitudes of the bodies. 
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2.1. Short-term numerical integrations 

Our first set of integrations was based on tlie complete equations of the or- 
bital motion of the satellites Remus and Romulus. These simulations we r e com - 
puted with the software NIMASTEP, presented in iDelsate fc Comperi (|201ll ). 
which is a home-made numerical software. It allows to integrate the osculating 
motion (using Cartesian coordinates) of an object considered as a point mass 
orbiting a homogeneous central body which rotates constantly around its princi- 
pal moment of inertia. It has been s u ccessfully tested and u s ed in many papers , 
as for e x ample iDelsate et al.l (IgOlOl). Ilemaitre et all (|2009( ). IValk et all (|2009l ). 



Delsate ( 2011 ). Compere et al.l ( 2012a ). To the aim of this work, the software 



has been improved in order to integrate the motion of two interacting satellites. 

The direct gravitational perturbations of the terrestrial and giant planets are 
negligible, as well as the orbital variations of Sylvia for the timespans considered 
here. We show in FigUthe importance of the two-body problem and of various 
perturbations on the acceleration of the satellite (only the radial component), 
in function of the relative distance from Sylvia, which is the ratio between 
the distance from the center of Sylvia and its equatorial radius. These results 
have been confirmed by numerical simulations. We consider in the following 
a Keplerian, eccentric and inclined orbit for Sylvia. The rotation of Sylvia is 
supposed to be constant and its axis of rotation, also constant, corresponds to 
its principal axis of inertia. 

We used the MEGNO indicator (jCincotta fc Simollioooh in order to distin- 
guish among regular and chaotic orbit. The MEGNO is a fast chaos indicator 
based on the numerical integration of a ta ngent vector, that has proved to be reli- 



able in a large class of dynamical systems (ICincotta et al.ll2003t iGozdziewski et al 
' af lBreiter et al.ll2005l:ICompere et al.l j2012a': 'Prouard et al."201lV 



2001L I20Q8S , , 

We numerically integrated the expressions oflCozdzicwski et al. (2001) along 
with the equations of motion. The initial tangent vector was chosen as random 
for each integrated orbit. Since the two satellites perturb each other, it is worth 
noting that the MEGNO indicator here was computed from the orbital evolution 
of the two satellites taken as a whole. The indicator thus represents the behavior 
of the complete system. 

We integrated a large set of orbits by varying the semi-major axis and the 
eccentricity of the two satellites, while keeping all the other initial variables as 
constant. The numerical integrations had a timespan of 20 years, until the satel- 
lites collided each other or with Sylvia, or in the case of ejection into heliocentric 
orbits. The results of the integrations are shown in Figl2] 

The actual position of Remus and Romulus lies in a very stable zone, even 
when considering their known incertitudes in semi-major axis and eccentricity 
(see Tabled]). The map clearly shows two types of instability zones (in black 
in the figure). One appears systematically for semi-major axis of Remus under 
400 km. Under this value, Remus quickly collides with Sylvia or is ejected from 
the system. The satellites have very low initial values of the eccentricity, which 
locates the orbits far from regions where overlapping of spin-orbit resonances is 
important. The width and shape of spin-orbit resonances in the general case of 



5 






* * * • 
















---GM 

— Sun 
Jupiter 

— Rad. Press, 
o Mars 


- 


* « « 




* « « « 
















* * . . 


* * * . 


* * * * . 


<^ 
* * + + 


<, 
* + + + 


'> <f 
* * * + + 


+ + + + 


+ * + + 




























































































o o o o 


o O o c 


o O O 


o o o o 


O O 


o o o o 


o o o o 


o o o o 


o o o o 


O O 


To o " ° 























1 2 3 4 5 6 7 8 9 10 11 12 



Relative distance from the center of the primary 



Figure 1: Order of magnitude of the accelerations of a satellite of Sylvia due to various 
perturbations, in function of the distance from Sylvia (expressed in radius of Sylvia). 



an orbiter have been investigated by Mvsen et al. ( 2006[) and Mvsen fc AksnesI 
A narrow vertical line of instability can thus be seen for Remus, at drem 
= 440 km, which corresponds to the spin-orbit resonance 3:1 (corresponding to 
the argument 3Arem — Aspm, where Xrem is the longitude of Remus and Xspin 
the spin angle of Sylvia). Higher-order resonances, like the 4:1 at arem = 530 
km are too weak to destabilize the orbit. 

The second type of instability zones corresponds to mean-motion resonances 
(MMRs) between the two satellites, at high semi-major axis in the Remus map, 
and low semi-major ax is in the Romulus map. The MMRs have a typical V- 
shape (see for example Gozdziewski et al. 2008lj: Bazso et al. 2010l ) which en- 
larges them for increasing eccentricity. These resonances can be seen clearlier 
in Figini where we modify the semi-major axis of the satellites in the range 
500 — 1500 km (Remus) and 1000 — 2000 km (Romulus). A same approach has 
been applied to investigate the dynamics of the giant planets by .Guzzo (2005L 
20061) . 



The MMRs correspond to the different lines, each one describing a resonance 
kitii ^ k2n2, where ni and n2 are the mean motions of the satellites, and fci, 
^2 integers. By neglecting the masses of the satellites with respect to the mass 
of Sylvia, the lines are described by the relation 02 — (fc2/fci)^^'^ai- The reso- 
nances overlap as the satellites have closer semi-major axis. The two important 
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Figure 2: MEGNO map of Remus (top) and Romulus (bottom) in semi-major axis and ec- 
centricity. Stable quasi-periodic orbits correspond to a color code of 2 (yellow), while chaotic 
orbits have a MEGNO > 2 (orange to black). A black cross indicates the actual position of 
the Remus-Romulus system. 



instability zones present on each side of the a2 — a-i hne represent initial or- 
bits leading to close encounters between the satellites, resulting in collisions or 
ejections from the system after a chaotic evolution. 

The actual position of the system is represented by a cross in the figure. The 
system thus lies in a very stable zone, bounded by the MMR 2:1 and 3:1, and 
shows a very regular evolution for at least 10'* years. Figure S] shows the period 
of Remus in function of the initial semi-major axis of Remus, while the initial 
semi-major axis of Romulus is equal to its actual value (thus investig ating a 
vertical line i n Figl3|). The period was obtained by a frequency analysis (jLaskar 
of the complex time series (a cos A, a sin A). 



Figure [5] shows the maximum eccentricity attained by the satellites over 20 
years, in function of their initial semi-major axes. The figure indicates that the 
eccentricity of Remus can possibly raise up to relatively large values inside the 
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Figure 3: MEGNO map in semi-major axis of both satellites. A black cross indicates the 
actual position of the Remus-Romulus system. 



major MMRs (up to 0.13 inside the MMR 2:1), and ten times more than the 
eccentricity of Romulus, which just reflects the fact that due to their respective 
masses, Remus is much more perturbed by Romulus than the reciprocal. 

We will now present the results obtained by investigating the secular evolu- 
tion of the system. 



2.2. Long-term numerical integrations 

The main limiting factor in the numerical integrations of the above model 
concerning CPU time is the fast rotation period of Sylvia (5.184 hours). We 
thus used equations of motion averaged over the spin angle of Sylvia, as well as 
over the mean anomalies of the satellites, to study the behavior of the system 
over longer timescales. 

We numerically integrated the Lagrange equations ( Brouwer fc Clemencd 



196lh where the averaged disturbing function (R) is the sum of different av- 



eraged perturbations (R) = (Rmut) + {Rou) + (Rq)- The mutual perturba- 
tion between the satellites (Rmut) is approximated by a fou rth-order expansion 
in eccentricity and inclination of their disturbing function ( Murrav fc Dermottl 
[2000). The oblateness perturbation (Rou) is represented by an expansion of 
fourth-order in eccentricity and inclinati on of the ob lateness disturbing func- 
tion, containing terms in J2, J|, and J4 (|yeras''2007 ). The additional se cular 
expression arising from C22 as derived bv lOe Saedeleer fc Henrardl ( 20061 ) was 
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Figure 4: Period of Remus obtained from frequency analysis, and location of the main mean- 
motion resonances. 



also taken into account, while proved to be of a very weak influence on the re- 
sults. Finall y, the solar pertu rbation (i?©) was modeled by using the analytical 
expansion of iBrumberg et al. which can be easily averaged over the mean 

anomalies of the satellite, while retaining the solar orbital evolution in spherical 
coordinates. The orbit of the Sun was modeled as a precessing orbit, where 
the three fundamental frequencies (tiq, tijQ, f^©) were taken from the AstDys 
database ( Knezevic fc Milani 2003() . 



The use of such averaged disturbing functions implies fixed semi-major axes 
for the satellites, as well as the suppression of the mean motion resonances from 
the dynamical system. The same orbits than in Figl3] (varying the semi-major 
axis of the satellites) are integrated over 6600 years. The results are shown in 
FiglH While the MEGNO could also be used in these integrations (by deriving 
the Lagrange equations over the orbital elements), we chose here to show the 
maximum eccentricity range attained by Romulus, which is also an indicator of 
the chaotic diffusion suffered by the satellite orbits. The white region indicates 
orbits for which the semi-major axis of Remus exceeds the one of Romulus, and 
are thus not numerically integrated. The color code range has been limited and 
chosen in order to magnify the different secular resonances. Vertical lines are 
visible and represent resonances belonging to the evection family, here between 
the pericenter frequency of Romulus and the mean motion of the Sun : 

kiw i2 /^2f^0, (2) 

with fci, ^2 integers. In particular, the so-called evection resonance (fci = k2 = 



9 



500 750 1000 1250 1500 

Initial Remus semi-major axis (km) 



750 1000 1250 1500 

Initial Remus semi-major axis (km) 



1200 1400 1600 1800 2000 

Initial Romulus semi-major axis (km) 




1200 1400 1600 1800 

Initial Romulus semi-major axis (km) 



Figure 5: Maximum eccentricity attained by Remus (left) and Romulus (right) over 20 years, 
in function of the initial semi-major axis of Remus and Romulus. The two important peaks 
seen in the two bottom panels correspond to the MMR 2:1 and 3:1. 



1) and its imp lications on th e dynamic s of satellites have been studied by 
many authors (IXouma fc Wisdo m 199 sl: iBreiteil bOOOt lYokovama et al.l bOOSi: 
Cuk &: Nesvornv 20101 : Frouard et al.ll2010() and is located in Figl6]at a Romu- 



lus semi-major axis of 1460 km. It is the most distinctive feature in the map. 
Inside the resonance, the orbit of Romulus becomes highly chaotic, its eccen- 
tricity being raised to 0.2. The other resonances from the evection family have 
a much weaker influence on the orbital evolutions. We stress that the evection 
resonances considered here are due to the oblateness of the asteroid, which ac- 
celerates considerably the precession period of w. The evection resonance can 
also be found much farther in semi-major axis, where the acceleration of the 
precession period is only due to the Sun. 

Despite the fact that the map shows the eccentricity range of Romulus, we 
can still see some effects due to the evection resonances for Remus. Such reso- 
nances are present in the upper right corner of the map; for example one starts 
at arom = 2000 km and a^em = 1250 km. These resonances should be approx- 
imately horizontal in the case where the mass of Remus is much greater than 
the one of Romulus, but here the perturbation of Romulus sensibly changes the 
frequency of the pericenter of Remus, with the effect of shaping these resonances 
as curves. 

Once again, the position of the actual system is shown on the map and lies 



10 



Romulus eccentricity range 




1000 1200 1400 1600 1800 2000 



Romulus semi-major axis (l<m) 

Figure 6; Maximum eccentricity range attained by Romulus over 6600 years, in function of 
the initial semi-major axis of the satellites. A green point indicates the actual position of the 
system. 



in a very stable zone. 



3. The effect of tides and BYORP 



We now investigate how tidal and BYORP effects can bring the system 
through the resona nces and what are t he tim esca les involved. 

As explained bv iGoldreich fc Saril (|2009l ) and lTavlor fc Margod (|2010l ). the 
tides raised by the primary on a secondary drive very efficiently the spin of the 
secondary Qg in synchronization with the mean motion n. We thus assume in 
the following that riRom = nnom and ^Rem = «i?.em- 

In a general way, the tidal equations describing t he variation of the semi- 
major axis and eccentrici ty of a satellite are, for e ^ 1 (jGoldreich fc Soteiill966 : 
Murrav fc Dermottllioooh : 



a — aign\\l^ — n)——^ I — -\ na, 



■T ■ (on a N57fc2p m, f Rp\ 

e = siqnililj, — 3n) — —-^ — ^ ne, 

^ ' % Qpmp \ a J 

for the tides raised on the planet by a (prograde) satellite, and 

5 



a = sign{!.ts — n)— \ — ) na, 



nis \ a 



(3) 
(4) 

(5) 
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• T _ _21 A|2s mp / R, 
2 Qs ms\ a 

for the tides raised on a prograde satellite by a planet, i? is the radius, m is 
the mass and the suffixes p and s represent respectively the primary body and 
the satellite. Note that the last expression in is non-null even if the satellite 
has its spin synchronized. The equations depend on k2 (which depends on the 
rigidity of the material) and on the dissipation function Q. 
Thus, when lOs = n, the tides raised on the planet by a satellite act to increase 
both the semi-major axis and the eccentricity, while those raised on the satellite 
by the planet are negligible for the semi-major axis and decrease the e ccentricity, 
some times completely counterbalancing the tides raised on the planet (jGoldreich 
Il963h . 



The tides depend on the parameters k2 and Q for each body, with 

i. 3/2 

1 + /i 

where p, is the dimensionless rigidity defined by : 

19fi 



M = 



2gpR' 



(7) 



(8) 



g is the gravity at the radius of the body {g — Gm/B?), p is its density, p is the 
rigidity of the material and G is the gravitational constant. Typical values of the 
rigidity are ^ 5x 10^°7V.w ~ ^ for a rocky body, or ^ 4x lO^A^.m"^ for an icy body 



( Murrav fc Dermottllioooh . iMarchis et aD (|2008al lbl) used 5 x 10^ N.m-^ consid- 
ering a moderately fractured aste roid for Sylvia, corresponding to fc2p ~ 0.015. 
However, iGoldreich ASarldlooi) showed how the effective rigidity of a rubble- 
pile can be considerably inferior to a monolith one, and gave the approximation 
P'Tubbie ^ [i^hvY^'^ where ey- is the yield strain and is taken to 10 ~^. This 



leads to oc R, although this was infirmed bv I Jacobson fc ScheeresI (2011) 



who found a relation k2p oc 1/R based on observations of the binary asteroid 
population. 

The value of the dissipation fonction Q is less known and is often chosen as 
~ 100 for monoliths (Goldreich & Sotter 1966). 

In this work we also take into acc ount the BYORP effect which arises from 
an as ymmetry of the satellite shapes ( Cuk fc Burnsll2005l: McMahon fc Scheeres 
2010h . We use the following equations from I Jacobson fc ScheeresI (|201ll ) : 
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^Q=10io, B=10-2 


AfQ=10", B=10-3 


Remus 


-47.8 Myr 


-45.7 Myr 


-453.2 Myr 


Romulus 


-200 Myr 


-182 Myr 


-1.8 Gyr 




1299.46 km 


1291.53 km 


1291.57 km 



Table 3: Dynamical ages of the satellites obtained from the tidal and BYORP evolution for 
different values of the tidal parameters and B. The satellites are considerated separately. 
The system is stopped when the satellites attain the semi-major axis limit a^rit = 400 km. 
The last row indicates the semi-major axis of Romulus when Remus attains a^rit- 



with ujii = -y/ AirGp/S and Hq — Fq / (og y 1 — Bq) where Fq is the solar radia- 
tion constant, a© and e© are the semi-major and eccentricity of the Sun. The 
BYORP effect depends on the parameter B which represents the deviation of 
the secondary with respect to a symmetric body. B is c ontained in the inter- 
val [0, 2] and is commonly taken to IQ-^. As explained in'jacobs on &: ScheeresI 
([2011;), the sign of and depend on the shape of the satellites. We assume 
first > and < 0. 

We can derive an approximate tidal age of the satellites by integrating nu- 
merically the tidal and BYORP equations backward in time. The satellites are 
considered separately. The evolution is stopped when the satellites attain the 
stability limit a = 400 km. We show the results for the representative tidal pa- 
rameters = [10", 10-3], [10^°, 10-3], [1010,10-2] in TableEl The tidal 
age of Romulus is of the same order as the one obtained by Taylor & Margot 
(2011) if we consider jiQ — 10^^. The semi-major axis of Romulus at the time 
when Remus attains the stability limit is also indicated, which reflects the tidal 
age of the actual system. 

The BYORP parameter B determines the amplitude of the BYORP effect, 
but has a reduced effect for the orbits considered here. We show in Figl7]the tidal 
evolution of the system with the three set of parameters. The evolution starts 
with the semi-major axis of Remus at 400 km up to the actual configuration, 
and is subsequently followed during 1 Gyr. The major resonances are shown 
in the figure. The evolution for the parameters fiQ=10^^, B=10~^ is merged 
with the one corresponding with /iQ=10i°, B=10~^, and ends shortly after the 
evection line. The satellites are on converging orbits (the orbit of Remus is 
expanding faster than the one of Romulus) for the three sets of parameters, and 
the system crosses the evection resonance before the mean- motion resonance 2:1 
in each case. 

As observed above, the BYORP effect is stronger for distant orbits. We can 
see in Figl7]that a larger value of B allows a faster evolution. It is interesting to 
note that it also reduces the distance between the two satellites, and thus drives 
the system towards the mean-motion resonances seen in Fig|31 and eventually 
towards highly chaotic zones. A smaller value of B keeps away the system from 
these zones as the integration time increases. A choice of i? — 10~^ also has 
the effect of decreasing the eccentricity of Remus after ~ 750 Myr. Its value 
monotically increases with the two other sets of parameters. 
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Figure 7: Paths followed by the system due to the effect of tides and BYORP with different 
values of the parameters /iQ and B. The actual position of the system is indicated by a 
cross. The solid line corresponds to the parameters fj,Q=10^^ , B=10~^ while the dashed line 
corresponds to fj,Q =10^'^ , B=10~^ . The evolution corresponding to the parameters fj,Q=lQ^^ , 
i?=10~^ is merged with the dashed line. See text for comments. 

Similarly, using inverse signs for and noticeably changes the evolutions 
only when B = 10~^. For this value, the orbits become divergent after 450 Myr, 
and the eccentricity of Remus can attain 0.08 after 1 Gyr, compared to 0.03 with 
our first choice of sign in the BYORP equations. The inverse signs have a very 
small effect on the dynamical ages of Table |31 

Finally, we can observe that the system will also encounter the evection 
resonance before the MMR 2:1 for the second choice of sign in the BYORP 
equation and the three set of parameters. We note that with B = 10^^, the 
system can even stay between the MMR 2:1 and 3:1 during 5 Gyr. 

4. Generic instability zones 

As we know now that the orbits of Remus and Romulus are located in a 
stable dynamical zone, it is interesting to know if it is also the case for the 
other known triple systems. We also extend this study to satellites of Main-belt 
asteroids, in the case where the masses of the satellites are significantly smaller 
than the mass of the asteroid. For greater relative masses of the satellites, one 
can no longer neglect the influence of their rotation and shape on their orbital 
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evolution (see Cuk &: Nesvornv 2010t ). 

The orbital elements of the asteroids have been taken from the AstDys 
database. Data concerning the physical parameters of the asteroids and as- 
sociated satellites, as well as the orbital elements of the satellites, have been 
taken from various sources; (22) Kalliope (Descamp s et alj 2008 : Marchis et al 



2008al). (45) Eugenia (jKaasalainen et al.|[2002 : Mar chis et al.l l2010l). (87) Svlvia 



dMarchis et a l 2005a). (93) Minerva ([Marchis et al.l201ll). (107) Camilla (iTorppa et al 



2003 : lMarchrs et al..2005cl). (121) Hermion e ([Marchis et al.l2005bl: iDescamps et al. 



2009 ), (130) Ele ktra (iMarchis et al.ll2008bl: DAM IT), (216) Kleopatra (iDescamps et al 



2011), (243) Ida (iPetit et al.lll997l). (283 ) Em ma (iMarchis et al.l2008bt DAMIT), 
(136108) Haumea (|Ragozzine fc Brownli2009l: iRabinowitz et al.ll2006l) . 

The determination of the location of the triple systems relatively to their 
mean motion resonances just requires the knowledge of the periods of the satel- 
lites. We show in Fig|S]the positions in semi-major axis of the five multiple sys- 
tems considered here : (45) Eugenia, (87) Sylvia, (93) Minerva, (216) Kleopatra, 
and (136108) Haumea. The semi-major axes are normalized by the Hill radius of 
the asteroid, computed with Rh = aast{'mast/{3MQ)y^^ , where Qast and mast 
are the semi-major axis and mass of the asteroid, and Mq is the solar mass. The 

Normalized positions of the multiple systems 



c 



X 



o 
'cc' 



03 

cn 

T3 
N 



0.018 
0.016 
0.014 
0.012 
0.01 
0.008 
0.006 
0.004 



1 1 


1 


1 ^ 






+ 


^ 1 1 


1 


1 



0.01 0.015 0.02 0.025 

Normalized semi-major axis (outer satellite) 



0.03 



Figure 8: Semi-major axis of the multiple systems normalized by the Hill radius of their 
asteroid. From left to right : (136108) Haumea, (216) Kleopatra, (93) Minerva, (87) Sylvia, 
(45) Eugenia. One of the possible tidal evolution of (87) Sylvia from Fig[7] is also shown 
(corresponding to the tidal parameters \i,Q = 10^", 5=10"'^). 



figure shows that all the multiple systems, except (216) Kleopatra, are between 
the mean motion resonances 2:1 and 3:1. Noting that the tidal and BYORP 
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evolutions move the systems outwards, it could be tempting to postulate that 
the MMR 2:1 prevents the evolution of the satellites from going beyond its lo- 
cation. In Figini we show the evolution of the eccentricity of Remus inside the 
MMR 2:1. The semi- major axes of the satellites are chosen in order to place 
their orbit inside the resonance according to the tidal evolution shown in Fig|5] 
(Orem=1067 km, arom=1694 km). The other initial elements are chosen equal 
to their actual values. The eccentricity of Romulus is much less perturbed and 
attains only 0.01 for the same timespan, and the inclinations of both satellites 
remain close to their initial values. Using the capture probability estimates of 
Dermott et al. I (ll988h . the satellites are assured to be captured in the resonance 



if their initial eccentricities are less than the critical values e*g„ — 0.094 and 
e* = 0.012. 
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Figure 9: Evolution of the eccentricity of Remus once tlie system is inside the mean motion 
resonance 2:1. 



The actual eccentricities of the satellites could thus indicate whether they 
suffered a chaotic evolution in the past. (216) Kleopatra is the only system 
that seems to have crossed the resonance but, unfortunately, we still have no 
informations about th e actual eccentricities of its satellites. 

[Ragozzi ne fc Brown (2009) studied the orbital and tidal evolution of (136108) 
Haumea, and argued that the passage of the system through the MMR 3:1 could 
explain the high value of eccentricity of its inner satellite Namaka (e = 0.249). 
However, it can be noted that the known eccentricites of the satellites in the 
three systems (45) Eugenia, (87) Sylvia and (93) Minerva are quite low; the 
maximum value is attained by Princess (the innermost satellite of Eugenia) 
with e = 0.069. This could indicate than the passage through the MMR 3:1 
is quite safe for the satellites. Although the location, dimensions and mass 
of Haumea make it a completely different object than the Main-Belt asteroids 
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studied here, with possibly different dynamical mecanisms, we can try to find 
another explanation for these differences in eccentricity. 

We show in the following the location of the satellites with respect to the 
evection resonance. In order to compute the semi-major axis corresponding to 
the evection a^vec for each satellite, we solve w = tiq for aet,ec- The frequency 
of the pericenter w is obtained with the corresponding Lagrange equation : 



1 



rVY~^d{R)_ 1 - cos(/) d{R) \ 



(11) 



where the averaged disturbing function (i?) — {Rou) + (^0) contains the solar 
in fluence and t he oblateness of the asteroid. (Rou) is taken from the expansion 
of Verad ( 2007 ) which keeps terms up to and sin^ /, and {Rq) is the classical 



Kozai' approximation. The eccentricity and inclination of the satellite are fixed 
values, but the solution Uevec still depends on the argument of the pericenter oj of 
the satellite, although we checked that this dependence is highly negligible. The 
oblateness coefficients J2 and J4 are determined from the shape of the ellipsoid 
approximation for each asteroid, using the formulation of iBovc 3 (ll997t) . 



An averaged disturbing function describing the perturbation of an additional 
satellite was used in the case of triple systems. The r esulting disturbing functio n 
(Rmut) was taken from the analytical expansion of Murray fc Dermott ( 200Cll ). 



In this case, the value of vu of each of the satellite was found while keeping fixed 
the semi-major of the other satellite. In the case of (87) Sylvia and (136108) 
Haumea, it was found that the frequency of pericenter of the innermost satellite 
could not be low enough to be commensurable with the mean motion of the 
Sun. In addition, the resonance can be found for the innermost satellite of 
(216) Kleopatra only when considering the first shape model, giving two values 
of Qevec ■ 915 km and 1456 km. However, we did not retain this system in 
the study, since the lack of information about the eccentricities of its satellites. 
Similarly, some binary systems have not been retained in this study due to the 
lack of information about the shape of their asteroid to this date. This is the 
case of (379) Huenna, (702) Alauda and (762) Pulkova. 

For some systems, both the shape of the asteroid from the ellipsoid ap- 
proximation, and a theorical value of J2 determined from the knowledge of the 
satellite orbit are known. The computation of aevec gives approximately similar 
results. 

The distance of the satellites from the evection resonance is shown in FigfTOl 
The uncertainties in semi-major axis and eccentricity correspond to the uncer- 
tainties on the knowledge of the orbits. The small number of systems available 
make statistical considerations rather ambiguous but we can try to make a few 
considerations, based on the assumption that the tidal and BYORP effects make 
the satellites evolve outwards. 

First we can consider the fact that satellites beyond the evection resonance 
are rather close to it. For a comparison, the Moon would be located at ^ 31 aevec 
in the figure (using for the Earth J2 = 1083 x lO^^ and J4 = -2 x 10"'^). Apart 
from observational biases, this could maybe point out the evection resonance 
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Figure 10: Semi-major axis and eccentricity of the satellites. The semi-major axis is normal- 
ized with the position of the evection resonance, represented by the vertical line. 



as a powerful way to eject satellites of asteroids. The effect of the evection 
resonanc e is already well-known for dramatically increasin g the eccentricity of a 
satellite (|Touma fc WisdomlHool ICuk fc Nesvornvll2010l ). As an example, we 
show in Figllll the highly chaotic evolution of Romulus inside the resonance. 
Its initial eccentricity, as well as the initial semi-major axis and eccentricity 
of Remus, have been given accordingly to the tidal and BYORP evolution of 
the system, with the parameters jiQ — 10^°, 5=10"'^. The larger number of 
satellites with very-low eccentricity located before the resonance could infer that 
hypothesis, which will be investigated in a following paper. 



5. Conclusion 

We studied the dynamics and stability of the (87) Sylvia triple system by 
using numerical integrations of the complete and averaged equations of motion. 
We used a shape of Sylvia derived from light-curves observations and up to 
C4,4, 5*4,4 for the complete integrations. The position of the actual system lies 
in a very stable zone. We showed the possible evolutions of the system trough 
tidal and BYORP effects, and showed that the system, currently lying between 
the mean- motion resonances 2:1 and 3:1, will likely evolve through the evection 
resonance before the MMR 2:1 in the future. The other known triple system 
considered here, except (216) Kleopatra, also lie between the MMR 2:1 and 3:1. 
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Figure 11: Eccentricity of Romulus inside tlie evection resonance. 

Finally, we show that the evection resonance could limit the outward evolution 
of the satellites. 
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